source('../env.R')
Using GitHub PAT from the git credential store.
Skipping install of 'clootl' from a github remote, the SHA1 (2ed1650b) has not changed since last install.
Use `force = TRUE` to force installation
community_data = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics_using_relative_abundance.csv'))
Rows: 308 Columns: 10── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (10): mntd_standard, mntd_actual, mass_fdiv_standard, mass_fdiv_actual, beak_width_fdiv_standard, beak_width_fdiv_actual, hwi_fdiv_...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(community_data)
colnames(community_data)
[1] "mntd_standard" "mntd_actual" "mass_fdiv_standard" "mass_fdiv_actual" "beak_width_fdiv_standard"
[6] "beak_width_fdiv_actual" "hwi_fdiv_standard" "hwi_fdiv_actual" "city_id" "urban_pool_size"
min(community_data$mntd_standard)
[1] -2.33692
max(community_data$mntd_standard)
[1] 2.328448
min(community_data$beak_width_fdiv_standard)
[1] -2.685152
max(community_data$beak_width_fdiv_standard)
[1] 1.931681
min(community_data$hwi_fdiv_standard)
[1] -2.200336
max(community_data$hwi_fdiv_standard)
[1] 2.333383
min(community_data$mass_fdiv_standard)
[1] -2.377212
max(community_data$mass_fdiv_standard)
[1] 2.1073
Join on realms
city_to_realm = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))
Rows: 337 Columns: 2── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
community_data_with_realm = left_join(community_data, city_to_realm)
Joining with `by = join_by(city_id)`
Cities as points
city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp'))) %>% left_join(community_data_with_realm)
Warning: st_centroid assumes attributes are constant over geometriesWarning: st_centroid does not give correct centroids for longitude/latitude dataJoining with `by = join_by(city_id)`
city_points_coords = st_coordinates(city_points)
city_points$latitude = city_points_coords[,1]
city_points$longitude = city_points_coords[,2]
world_map = read_country_boundaries()
Load community data, and create long format version
communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))
Rows: 2428 Columns: 7── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, ebird_species_name, seasonal, presence, origin
dbl (2): city_id, relative_abundance_proxy
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
communities
community_summary = communities %>% group_by(city_id) %>% summarise(regional_pool_size = n(), urban_pool_size = sum(relative_abundance_proxy > 0))
community_summary
Load trait data
traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_ebird.csv'))
Rows: 332 Columns: 10── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): ebird_species_name, habitat, trophic_level, trophic_niche, primary_lifestyle
dbl (5): beak_width, hwi, mass, habitat_density, migration
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(traits)
Load realm geo
resolve = read_resolve()
head(resolve)
Simple feature collection with 6 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -162.1547 ymin: -69.55876 xmax: 158.6167 ymax: 61.53428
Geodetic CRS: WGS 84
test_required_values = function(name, df) {
cat(paste(
test_value_wilcox(paste(name, 'MNTD'), df$mntd_standard),
test_value_wilcox(paste(name, 'Beak Gape FDiv'), df$beak_width_fdiv_standard),
test_value_wilcox(paste(name, 'HWI FDiv'), df$hwi_fdiv_standard),
test_value_wilcox(paste(name, 'Mass FDiv'), df$mass_fdiv_standard),
paste('N', nrow(df)),
sep = "\n"))
}
test_required_values('Global', community_data_with_realm)
Global MNTD median -0.36 ***
Global Beak Gape FDiv median 0.02
Global HWI FDiv median 0.39 **
Global Mass FDiv median 0.29 ***
N 308
unique(community_data_with_realm$core_realm)
[1] "Nearctic" "Neotropic" "Palearctic" "Afrotropic" "Indomalayan" "Australasia"
test_required_values('Nearctic', community_data_with_realm[community_data_with_realm$core_realm == 'Nearctic',])
Nearctic MNTD median 0.67 *
Nearctic Beak Gape FDiv median 0.29
Nearctic HWI FDiv median -0.8 ***
Nearctic Mass FDiv median -0.26
N 46
test_required_values('Neotropic', community_data_with_realm[community_data_with_realm$core_realm == 'Neotropic',])
Neotropic MNTD median 0.03
Neotropic Beak Gape FDiv median -0.44 ***
Neotropic HWI FDiv median -0.31
Neotropic Mass FDiv median 0.33 *
N 64
test_required_values('Palearctic', community_data_with_realm[community_data_with_realm$core_realm == 'Palearctic',])
Palearctic MNTD median 0.13
Palearctic Beak Gape FDiv median 1.25 ***
Palearctic HWI FDiv median -0.39
Palearctic Mass FDiv median 0.01
N 72
test_required_values('Afrotropic', community_data_with_realm[community_data_with_realm$core_realm == 'Afrotropic',])
Afrotropic MNTD median -1.28 *
Afrotropic Beak Gape FDiv median -0.56
Afrotropic HWI FDiv median 0.15
Afrotropic Mass FDiv median -0.95
N 9
test_required_values('Indomalayan', community_data_with_realm[community_data_with_realm$core_realm == 'Indomalayan',])
Indomalayan MNTD median -0.64 ***
Indomalayan Beak Gape FDiv median -0.68 ***
Indomalayan HWI FDiv median 1.11 ***
Indomalayan Mass FDiv median 0.83 ***
N 111
test_required_values('Australasia', community_data_with_realm[community_data_with_realm$core_realm == 'Australasia',])
Australasia MNTD median -1.39
Australasia Beak Gape FDiv median -0.75
Australasia HWI FDiv median 0.77
Australasia Mass FDiv median -0.96
N 6
communities %>%
left_join(city_to_realm) %>%
mutate(family = gsub( "_.*$", "", ebird_species_name)) %>%
dplyr::select(family, core_realm) %>%
distinct() %>%
arrange(core_realm)
Joining with `by = join_by(city_id)`
communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))
Rows: 2428 Columns: 7── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, ebird_species_name, seasonal, presence, origin
dbl (2): city_id, relative_abundance_proxy
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
city_introduced_species = communities %>% group_by(city_id) %>% summarise(number_of_species = n()) %>% left_join(
communities %>% group_by(city_id) %>% filter(origin == 'Introduced') %>% summarise(number_of_introduced_species = n())
) %>% replace_na(list(number_of_introduced_species = 0))
Joining with `by = join_by(city_id)`
community_data_with_introductions = left_join(community_data, city_introduced_species)
Joining with `by = join_by(city_id)`
community_data_with_introductions$has_introduced_species = community_data_with_introductions$number_of_introduced_species > 0
community_data_with_introductions
community_data_with_introductions[,c('mntd_standard', 'has_introduced_species')]
community_data_with_introductions %>% group_by(has_introduced_species) %>% summarise(
total_cities = n(),
mean_mntd_std = mean(mntd_standard, na.rm = T),
median_mntd_std = median(mntd_standard, na.rm = T),
sd_mntd_std = sd(mntd_standard, na.rm = T),
mean_mass_fdiv_std = mean(mass_fdiv_standard, na.rm = T),
median_mass_fdiv_std = median(mass_fdiv_standard, na.rm = T),
sd_mass_fdiv_std = sd(mass_fdiv_standard, na.rm = T),
mean_gape_width_fdiv_std = mean(beak_width_fdiv_standard, na.rm = T),
median_gape_width_fdiv_std = median(beak_width_fdiv_standard, na.rm = T),
sd_gape_width_fdiv_std = sd(beak_width_fdiv_standard, na.rm = T),
mean_handwing_index_fdiv_std = mean(hwi_fdiv_standard, na.rm = T),
median_handwing_index_fdiv_std = median(hwi_fdiv_standard, na.rm = T),
sd_handwing_index_fdiv_std = sd(hwi_fdiv_standard, na.rm = T)
)
ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = mntd_standard)) + geom_boxplot()
wilcox.test(mntd_standard ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')
Wilcoxon rank sum test with continuity correction
data: mntd_standard by has_introduced_species
W = 7925, p-value = 0.00001285
alternative hypothesis: true location shift is not equal to 0
There is a significant difference between the response of cities with introduced species (0.53±0.27) and those without (0.47±0.19) (p-value = 0.02).
ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = mass_fdiv_standard)) + geom_boxplot()
wilcox.test(mass_fdiv_standard ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')
Wilcoxon rank sum test with continuity correction
data: mass_fdiv_standard by has_introduced_species
W = 15028, p-value = 0.0000006706
alternative hypothesis: true location shift is not equal to 0
There is a significant difference between the response of cities with introduced species (0.57±0.27) and those without (0.73±0.24) (p < 0.0001)
ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = beak_width_fdiv_standard)) + geom_boxplot()
wilcox.test(beak_width_fdiv_standard ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')
Wilcoxon rank sum test with continuity correction
data: beak_width_fdiv_standard by has_introduced_species
W = 8662, p-value = 0.0006884
alternative hypothesis: true location shift is not equal to 0
There is NOT a significant difference between the response of cities with introduced species (0.61±0.30) and those without (0.56±0.27)
ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = hwi_fdiv_standard)) + geom_boxplot()
wilcox.test(hwi_fdiv_standard ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')
Wilcoxon rank sum test with continuity correction
data: hwi_fdiv_standard by has_introduced_species
W = 17606, p-value < 0.00000000000000022
alternative hypothesis: true location shift is not equal to 0
There is a significant difference between the response of cities with introduced species (0.49±0.30) and those without (0.79±0.21) (p < 0.0001)
geography = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'geography.csv'))
Rows: 342 Columns: 26── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (26): city_id, city_avg_ndvi, city_avg_elevation, city_avg_temp, city_avg_min_monthly_temp, city_avg_max_monthly_temp, city_avg_mon...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
names(geography)
[1] "city_id" "city_avg_ndvi" "city_avg_elevation" "city_avg_temp"
[5] "city_avg_min_monthly_temp" "city_avg_max_monthly_temp" "city_avg_monthly_temp" "city_avg_rainfall"
[9] "city_avg_max_monthly_rainfall" "city_avg_min_monthly_rainfall" "city_avg_soil_moisture" "city_max_elev"
[13] "city_min_elev" "city_elev_range" "region_20km_avg_ndvi" "region_20km_avg_elevation"
[17] "region_20km_avg_soil_moisture" "region_20km_max_elev" "region_20km_min_elev" "region_20km_elev_range"
[21] "region_50km_avg_ndvi" "region_50km_avg_elevation" "region_50km_avg_soil_moisture" "region_50km_max_elev"
[25] "region_50km_min_elev" "region_50km_elev_range"
analysis_data = community_data_with_realm[,c('city_id', 'mntd_standard', 'mass_fdiv_standard', 'beak_width_fdiv_standard', 'hwi_fdiv_standard', 'core_realm')] %>%
left_join(city_points[,c('city_id', 'latitude', 'longitude')]) %>%
left_join(community_data_with_introductions[,c('city_id', 'has_introduced_species')]) %>%
left_join(geography)
Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`
analysis_data$abs_latitude = abs(analysis_data$latitude)
analysis_data$core_realm = factor(analysis_data$core_realm, levels = c('Palearctic', 'Nearctic', 'Neotropic', 'Afrotropic', 'Indomalayan', 'Australasia', 'Oceania'))
analysis_data$has_introduced_species = factor(analysis_data$has_introduced_species, level = c('TRUE', 'FALSE'), labels = c('Introduced species', 'No introduced species'))
model_data = function(df, dependant_var) {
df[,c(dependant_var, 'core_realm', 'abs_latitude', 'longitude', 'has_introduced_species', 'city_avg_ndvi', 'city_avg_elevation', 'city_avg_temp', 'city_avg_min_monthly_temp', 'city_avg_max_monthly_temp', 'city_avg_monthly_temp', 'city_avg_rainfall', 'city_avg_max_monthly_rainfall', 'city_avg_min_monthly_rainfall', 'city_avg_soil_moisture', 'city_max_elev', 'city_min_elev', 'city_elev_range', 'region_20km_avg_ndvi', 'region_20km_avg_elevation', 'region_20km_avg_soil_moisture', 'region_20km_max_elev', 'region_20km_min_elev', 'region_20km_elev_range', 'region_50km_avg_ndvi', 'region_50km_avg_elevation', 'region_50km_avg_soil_moisture', 'region_50km_max_elev', 'region_50km_min_elev', 'region_50km_elev_range')]
}
model_data(analysis_data, 'mntd_standard')
names(analysis_data)
[1] "city_id" "mntd_standard" "mass_fdiv_standard" "beak_width_fdiv_standard"
[5] "hwi_fdiv_standard" "core_realm" "latitude" "longitude"
[9] "geometry" "has_introduced_species" "city_avg_ndvi" "city_avg_elevation"
[13] "city_avg_temp" "city_avg_min_monthly_temp" "city_avg_max_monthly_temp" "city_avg_monthly_temp"
[17] "city_avg_rainfall" "city_avg_max_monthly_rainfall" "city_avg_min_monthly_rainfall" "city_avg_soil_moisture"
[21] "city_max_elev" "city_min_elev" "city_elev_range" "region_20km_avg_ndvi"
[25] "region_20km_avg_elevation" "region_20km_avg_soil_moisture" "region_20km_max_elev" "region_20km_min_elev"
[29] "region_20km_elev_range" "region_50km_avg_ndvi" "region_50km_avg_elevation" "region_50km_avg_soil_moisture"
[33] "region_50km_max_elev" "region_50km_min_elev" "region_50km_elev_range" "abs_latitude"
all_explanatories = c(
'abs_latitude', 'latitude', 'longitude',
'has_introduced_species',
'city_avg_ndvi', 'city_avg_elevation', 'city_avg_temp', 'city_avg_min_monthly_temp', 'city_avg_max_monthly_temp',
'city_avg_monthly_temp', 'city_avg_rainfall', 'city_avg_max_monthly_rainfall', 'city_avg_min_monthly_rainfall',
'city_avg_soil_moisture', 'city_max_elev', 'city_min_elev', 'city_elev_range',
'region_20km_avg_ndvi', 'region_20km_avg_elevation', 'region_20km_avg_soil_moisture', 'region_20km_max_elev',
'region_20km_min_elev', 'region_20km_elev_range',
'region_50km_avg_ndvi', 'region_50km_avg_elevation', 'region_50km_avg_soil_moisture', 'region_50km_max_elev',
'region_50km_min_elev', 'region_50km_elev_range',
'core_realmAfrotropic', 'core_realmAustralasia', 'core_realmIndomalayan', 'core_realmNearctic', 'core_realmNeotropic', 'core_realmPalearctic')
type_labels = function(p) {
explanatory_levels = all_explanatories[all_explanatories %in% p$explanatory]
p$explanatory <- factor(p$explanatory, levels = explanatory_levels)
p$type <- 'Realm'
p$type[p$explanatory %in% c('city_avg_ndvi', 'city_avg_elevation', 'city_avg_temp', 'city_avg_min_monthly_temp', 'city_avg_max_monthly_temp',
'city_avg_monthly_temp', 'city_avg_rainfall', 'city_avg_max_monthly_rainfall', 'city_avg_min_monthly_rainfall',
'city_avg_soil_moisture', 'city_max_elev', 'city_min_elev', 'city_elev_range')] <- 'City geography'
p$type[p$explanatory %in% c('region_50km_avg_ndvi', 'region_50km_avg_elevation', 'region_50km_avg_soil_moisture', 'region_50km_max_elev',
'region_50km_min_elev', 'region_50km_elev_range')] <- 'Regional (50km) geography'
p$type[p$explanatory %in% c('region_20km_avg_ndvi', 'region_20km_avg_elevation', 'region_20km_avg_soil_moisture', 'region_20km_max_elev',
'region_20km_min_elev', 'region_20km_elev_range')] <- 'Regional (20km) geography'
p$type[p$explanatory %in% c('abs_latitude', 'longitude')] <- 'Spatial'
p
}
explanatory_labels = c(
'has_introduced_species'='Has introduced species',
'city_avg_ndvi'='Average NDVI',
'city_avg_elevation'='Average elevation',
'city_avg_temp'='Average temperature',
'city_avg_min_monthly_temp'='Average minimum monthly temperature',
'city_avg_max_monthly_temp'='Average maximum monthly temperature',
'city_avg_monthly_temp'='Average monthly temperature',
'city_avg_rainfall'='Average rainfall',
'city_avg_max_monthly_rainfall'='Average maximum monthly rainfall',
'city_avg_min_monthly_rainfall'='Average minimum monthly rainfall',
'city_avg_soil_moisture'='Average soil moisture',
'city_max_elev'='Maximum elevation',
'city_min_elev'='Minimum elevation',
'city_elev_range'='Elevation range',
'region_20km_avg_ndvi'='Average NDVI',
'region_20km_avg_elevation'='Average elevation',
'region_20km_avg_soil_moisture'='Average soil moisture',
'region_20km_max_elev'='Maximum elevation',
'region_20km_min_elev'='Minimum elevation',
'region_20km_elev_range'='Elevation range',
'region_50km_avg_ndvi'='Average NDVI',
'region_50km_avg_elevation'='Average elevation',
'region_50km_avg_soil_moisture'='Average soil moisture',
'region_50km_max_elev'='Maximum elevation',
'region_50km_min_elev'='Minimum elevation',
'region_50km_elev_range'='Elevation range',
'abs_latitude' = 'Absolute latitude',
'latitude' = 'Latitude',
'longitude' = 'Longitude',
'core_realmAfrotropic' = 'Afrotropical',
'core_realmAustralasia' = 'Austaliasian',
'core_realmIndomalayan' = 'Indomalayan',
'core_realmNearctic' = 'Nearctic',
'core_realmNeotropic' = 'Neotropical',
'core_realmPalearctic' = 'Palearctic',
'core_realmOceania' = 'Oceanical')
geom_map = function(map_sf, title) {
norm_mntd_analysis_geo = ggplot() +
geom_sf(data = world_map, aes(geometry = geometry)) +
map_sf +
standardised_colours_scale +
labs(colour = 'Standardised\nResponse') +
theme_bw() +
theme(legend.position="bottom")
}
# Taken from MuMIN package
# https://rdrr.io/cran/MuMIn/src/R/averaging.R
# https://rdrr.io/cran/MuMIn/src/R/model.avg.R
.coefarr.avg <-
function(cfarr, weight, revised.var, full, alpha) {
weight <- weight / sum(weight)
nCoef <- dim(cfarr)[3L]
if(full) {
nas <- is.na(cfarr[, 1L, ]) & is.na(cfarr[, 2L, ])
cfarr[, 1L, ][nas] <- cfarr[, 2L, ][nas] <- 0
#cfarr[, 1L:2L, ][is.na(cfarr[, 1L:2L, ])] <- 0
if(!all(is.na(cfarr[, 3L, ])))
cfarr[ ,3L, ][is.na(cfarr[ , 3L, ])] <- Inf
}
avgcoef <- array(dim = c(nCoef, 5L),
dimnames = list(dimnames(cfarr)[[3L]], c("Estimate",
"Std. Error", "Adjusted SE", "Lower CI", "Upper CI")))
for(i in seq_len(nCoef))
avgcoef[i, ] <- par.avg(cfarr[, 1L, i], cfarr[, 2L, i], weight,
df = cfarr[, 3L, i], alpha = alpha, revised.var = revised.var)
avgcoef[is.nan(avgcoef)] <- NA
return(avgcoef)
}
.makecoefmat <- function(cf) {
no.ase <- all(is.na(cf[, 3L]))
z <- abs(cf[, 1L] / cf[, if(no.ase) 2L else 3L])
pval <- 2 * pnorm(z, lower.tail = FALSE)
cbind(cf[, if(no.ase) 1L:2L else 1L:3L, drop = FALSE],
`z value` = z, `Pr(>|z|)` = zapsmall(pval))
}
# Generate model selections using lmer, dredge, and model.avg
# `forumla` : a two-sided linear formula object describing both the fixed-effects and random-effects part of the model
# `data` : the data frame containing the variables from the formula
# `aic_delta` : the AIC delta to use for selecting models in model average
model_average <- function(formula, data, aic_delta = 20) {
model <- lm(
formula,
data=data
)
dredge_result <- dredge(model)
summary(model.avg(dredge_result, subset = delta < aic_delta))
}
# Create a summary data frame containing the selected variables from a model
# `model_sum` : The model summary output from `model_average`
model_summary <- function(model_sum) {
.column_name <- function(postfix) {
postfix
}
# just return the estimate and p value
weight <- model_sum$msTable[, 5L]
coefmat.full <- as.data.frame(.makecoefmat(.coefarr.avg(model_sum$coefArray, weight,
attr(model_sum, "revised.var"), TRUE, 0.05)))
coefmat.subset <-
as.data.frame(.makecoefmat(.coefarr.avg(model_sum$coefArray, weight,
attr(model_sum, "revised.var"), FALSE, 0.05)))
coefmat.subset <- coefmat.subset[-c(1), c(1, 2, 5)]
names(coefmat.subset) <- c(.column_name("estimate"), .column_name("error"), .column_name("p"))
coefmat.subset <- tibble::rownames_to_column(coefmat.subset, "explanatory")
coefmat.subset$model = 'subset'
coefmat.full <- coefmat.full[-c(1), c(1, 2, 5)]
names(coefmat.full) <- c(.column_name("estimate"), .column_name("error"), .column_name("p"))
coefmat.full <- tibble::rownames_to_column(coefmat.full, "explanatory")
coefmat.full$model = 'full'
rbind(coefmat.full, coefmat.subset)
}
formula_from_vsurp = function(predictors, dependent, vsurp_result) {
as.formula(paste(dependent, paste(names(predictors[,vsurp_result$varselect.interp]), collapse="+"), sep = "~"))
}
plot_vsurp_result = function(result_table) {
p = result_table[result_table$model == 'full',]
p = type_labels(p)
ggplot(p, aes(y=explanatory, x=estimate, colour = type)) +
geom_line() +
geom_point()+
geom_errorbar(aes(xmin=estimate-error, xmax=estimate+error), width=.2,
position=position_dodge(0.05)) +
scale_y_discrete(
limits = rev(levels(p$explanatory)),
labels = explanatory_labels) +
scale_colour_manual(values = c(colour1, colour2, colour3, colour4, colour5, colour6), breaks = c('Realm', 'City geography', 'Regional (50km) geography', 'Regional (20km) geography', 'Spatial')) +
theme_bw() +
geom_vline(xintercept=0, linetype="dotted") +
guides(colour=guide_legend(title="Predictor type")) + xlab('Difference in response from 0\nhabitat filtering (< 0) and competitive interactions (> 0)\n± Standard Error') + ylab('Predictor') +
theme(legend.justification = "top")
}
std_mntd_analysis_geo_plot = geom_map(geom_sf(data = analysis_data, aes(color = mntd_standard, geometry = geometry)), 'MNTD')
std_mntd_analysis_geo_plot
std_mntd_analysis_data = model_data(analysis_data[!is.na(analysis_data$mntd_standard),], 'mntd_standard')
std_mntd_analysis_predictors = std_mntd_analysis_data[,-1]
std_mntd_analysis_interp = VSURF(x = std_mntd_analysis_predictors, y = std_mntd_analysis_data$mntd_standard)
Thresholding step
Estimated computational time (on one core): 14 sec.
|
| | 0%
|
|====== | 5%
|
|============= | 10%
|
|=================== | 15%
|
|========================== | 20%
|
|================================ | 25%
|
|======================================= | 30%
|
|============================================= | 35%
|
|==================================================== | 40%
|
|========================================================== | 45%
|
|================================================================ | 50%
|
|======================================================================= | 55%
|
|============================================================================= | 60%
|
|==================================================================================== | 65%
|
|========================================================================================== | 70%
|
|================================================================================================= | 75%
|
|======================================================================================================= | 80%
|
|============================================================================================================== | 85%
|
|==================================================================================================================== | 90%
|
|=========================================================================================================================== | 95%
|
|=================================================================================================================================| 100%
Interpretation step (on 29 variables)
Estimated computational time (on one core): between 1.5 sec. and 28.4 sec.
|
| | 0%
|
|==== | 3%
|
|========= | 7%
|
|============= | 10%
|
|================== | 14%
|
|====================== | 17%
|
|=========================== | 21%
|
|=============================== | 24%
|
|==================================== | 28%
|
|======================================== | 31%
|
|============================================ | 34%
|
|================================================= | 38%
|
|===================================================== | 41%
|
|========================================================== | 45%
|
|============================================================== | 48%
|
|=================================================================== | 52%
|
|======================================================================= | 55%
|
|============================================================================ | 59%
|
|================================================================================ | 62%
|
|===================================================================================== | 66%
|
|========================================================================================= | 69%
|
|============================================================================================= | 72%
|
|================================================================================================== | 76%
|
|====================================================================================================== | 79%
|
|=========================================================================================================== | 83%
|
|=============================================================================================================== | 86%
|
|==================================================================================================================== | 90%
|
|======================================================================================================================== | 93%
|
|============================================================================================================================= | 97%
|
|=================================================================================================================================| 100%
Prediction step (on 11 variables)
Maximum estimated computational time (on one core): 4.2 sec.
|
| | 0%
|
|============ | 9%
|
|======================= | 18%
|
|=================================== | 27%
|
|=============================================== | 36%
|
|=========================================================== | 45%
|
|====================================================================== | 55%
|
|================================================================================== | 64%
|
|============================================================================================== | 73%
|
|========================================================================================================== | 82%
|
|===================================================================================================================== | 91%
|
|=================================================================================================================================| 100%
names(std_mntd_analysis_predictors[,std_mntd_analysis_interp$varselect.interp])
[1] "core_realm" "longitude" "city_avg_max_monthly_temp" "city_avg_temp"
[5] "abs_latitude" "city_avg_max_monthly_rainfall" "city_avg_monthly_temp" "city_avg_min_monthly_temp"
[9] "city_avg_ndvi" "city_avg_rainfall" "city_max_elev"
std_mntd_analysis_formula = formula_from_vsurp(std_mntd_analysis_predictors, "mntd_standard", std_mntd_analysis_interp)
std_mntd_analysis_result <- model_average(std_mntd_analysis_formula, std_mntd_analysis_data)
Fixed term is "(Intercept)"
std_mntd_analysis_result_table = model_summary(std_mntd_analysis_result)
std_mntd_analysis_result_table
std_mntd_analysis_pred_plot = plot_vsurp_result(std_mntd_analysis_result_table)
std_mntd_analysis_pred_plot
std_gape_fdiv_analysis_geo_plot = geom_map(geom_sf(data = analysis_data, aes(color = beak_width_fdiv_standard, geometry = geometry)), 'Beak Width FDiv')
std_gape_fdiv_analysis_geo_plot
std_gape_fdiv_analysis_data = model_data(analysis_data[!is.na(analysis_data$beak_width_fdiv_standard),], 'beak_width_fdiv_standard')
std_gape_fdiv_analysis_predictors = std_gape_fdiv_analysis_data[,-1]
std_gape_fdiv_analysis_interp = VSURF(x = std_gape_fdiv_analysis_predictors, y = std_gape_fdiv_analysis_data$beak_width_fdiv_standard)
Thresholding step
Estimated computational time (on one core): 13.3 sec.
|
| | 0%
|
|====== | 5%
|
|============= | 10%
|
|=================== | 15%
|
|========================== | 20%
|
|================================ | 25%
|
|======================================= | 30%
|
|============================================= | 35%
|
|==================================================== | 40%
|
|========================================================== | 45%
|
|================================================================ | 50%
|
|======================================================================= | 55%
|
|============================================================================= | 60%
|
|==================================================================================== | 65%
|
|========================================================================================== | 70%
|
|================================================================================================= | 75%
|
|======================================================================================================= | 80%
|
|============================================================================================================== | 85%
|
|==================================================================================================================== | 90%
|
|=========================================================================================================================== | 95%
|
|=================================================================================================================================| 100%
Interpretation step (on 29 variables)
Estimated computational time (on one core): between 5.5 sec. and 25.8 sec.
|
| | 0%
|
|==== | 3%
|
|========= | 7%
|
|============= | 10%
|
|================== | 14%
|
|====================== | 17%
|
|=========================== | 21%
|
|=============================== | 24%
|
|==================================== | 28%
|
|======================================== | 31%
|
|============================================ | 34%
|
|================================================= | 38%
|
|===================================================== | 41%
|
|========================================================== | 45%
|
|============================================================== | 48%
|
|=================================================================== | 52%
|
|======================================================================= | 55%
|
|============================================================================ | 59%
|
|================================================================================ | 62%
|
|===================================================================================== | 66%
|
|========================================================================================= | 69%
|
|============================================================================================= | 72%
|
|================================================================================================== | 76%
|
|====================================================================================================== | 79%
|
|=========================================================================================================== | 83%
|
|=============================================================================================================== | 86%
|
|==================================================================================================================== | 90%
|
|======================================================================================================================== | 93%
|
|============================================================================================================================= | 97%
|
|=================================================================================================================================| 100%
Prediction step (on 11 variables)
Maximum estimated computational time (on one core): 4.2 sec.
|
| | 0%
|
|============ | 9%
|
|======================= | 18%
|
|=================================== | 27%
|
|=============================================== | 36%
|
|=========================================================== | 45%
|
|====================================================================== | 55%
|
|================================================================================== | 64%
|
|============================================================================================== | 73%
|
|========================================================================================================== | 82%
|
|===================================================================================================================== | 91%
|
|=================================================================================================================================| 100%
names(std_gape_fdiv_analysis_predictors[,std_gape_fdiv_analysis_interp$varselect.interp])
[1] "abs_latitude" "core_realm" "longitude" "city_avg_max_monthly_temp"
[5] "city_avg_min_monthly_temp" "city_avg_temp" "city_avg_max_monthly_rainfall" "city_avg_monthly_temp"
[9] "city_avg_min_monthly_rainfall" "city_max_elev" "region_50km_max_elev"
std_gape_fdiv_analysis_formula = formula_from_vsurp(std_gape_fdiv_analysis_predictors, "beak_width_fdiv_standard", std_gape_fdiv_analysis_interp)
std_gape_fdiv_analysis_result <- model_average(std_gape_fdiv_analysis_formula, std_gape_fdiv_analysis_data)
Fixed term is "(Intercept)"
std_gape_fdiv_analysis_result_table = model_summary(std_gape_fdiv_analysis_result)
std_gape_fdiv_analysis_result_table
std_gape_fdiv_analysis_pred_plot = plot_vsurp_result(std_gape_fdiv_analysis_result_table)
std_gape_fdiv_analysis_pred_plot
std_hwi_fdiv_analysis_geo_plot = geom_map(geom_sf(data = analysis_data, aes(color = hwi_fdiv_standard, geometry = geometry)), 'HWI FDiv')
std_hwi_fdiv_analysis_geo_plot
std_hwi_fdiv_analysis_data = model_data(analysis_data[!is.na(analysis_data$hwi_fdiv_standard),], 'hwi_fdiv_standard')
std_hwi_fdiv_analysis_predictors = std_hwi_fdiv_analysis_data[,-1]
std_hwi_fdiv_analysis_interp = VSURF(x = std_hwi_fdiv_analysis_predictors, y = std_hwi_fdiv_analysis_data$hwi_fdiv_standard)
Thresholding step
Estimated computational time (on one core): 13.1 sec.
|
| | 0%
|
|====== | 5%
|
|============= | 10%
|
|=================== | 15%
|
|========================== | 20%
|
|================================ | 25%
|
|======================================= | 30%
|
|============================================= | 35%
|
|==================================================== | 40%
|
|========================================================== | 45%
|
|================================================================ | 50%
|
|======================================================================= | 55%
|
|============================================================================= | 60%
|
|==================================================================================== | 65%
|
|========================================================================================== | 70%
|
|================================================================================================= | 75%
|
|======================================================================================================= | 80%
|
|============================================================================================================== | 85%
|
|==================================================================================================================== | 90%
|
|=========================================================================================================================== | 95%
|
|=================================================================================================================================| 100%
Interpretation step (on 29 variables)
Estimated computational time (on one core): between 1.5 sec. and 26.1 sec.
|
| | 0%
|
|==== | 3%
|
|========= | 7%
|
|============= | 10%
|
|================== | 14%
|
|====================== | 17%
|
|=========================== | 21%
|
|=============================== | 24%
|
|==================================== | 28%
|
|======================================== | 31%
|
|============================================ | 34%
|
|================================================= | 38%
|
|===================================================== | 41%
|
|========================================================== | 45%
|
|============================================================== | 48%
|
|=================================================================== | 52%
|
|======================================================================= | 55%
|
|============================================================================ | 59%
|
|================================================================================ | 62%
|
|===================================================================================== | 66%
|
|========================================================================================= | 69%
|
|============================================================================================= | 72%
|
|================================================================================================== | 76%
|
|====================================================================================================== | 79%
|
|=========================================================================================================== | 83%
|
|=============================================================================================================== | 86%
|
|==================================================================================================================== | 90%
|
|======================================================================================================================== | 93%
|
|============================================================================================================================= | 97%
|
|=================================================================================================================================| 100%
Prediction step (on 6 variables)
Maximum estimated computational time (on one core): 1.8 sec.
|
| | 0%
|
|====================== | 17%
|
|=========================================== | 33%
|
|================================================================ | 50%
|
|====================================================================================== | 67%
|
|============================================================================================================ | 83%
|
|=================================================================================================================================| 100%
names(std_hwi_fdiv_analysis_predictors[,std_hwi_fdiv_analysis_interp$varselect.interp])
[1] "core_realm" "city_avg_temp" "city_avg_max_monthly_temp" "region_50km_min_elev"
[5] "longitude" "abs_latitude"
std_hwi_fdiv_analysis_formula = formula_from_vsurp(std_hwi_fdiv_analysis_predictors, "hwi_fdiv_standard", std_hwi_fdiv_analysis_interp)
std_hwi_fdiv_analysis_result <- model_average(std_hwi_fdiv_analysis_formula, std_hwi_fdiv_analysis_data)
Fixed term is "(Intercept)"
std_hwi_fdiv_analysis_result_table = model_summary(std_hwi_fdiv_analysis_result)
std_hwi_fdiv_analysis_result_table
std_hwi_fdiv_analysis_pred_plot = plot_vsurp_result(std_hwi_fdiv_analysis_result_table)
std_hwi_fdiv_analysis_pred_plot
std_mass_fdiv_analysis_geo_plot = geom_map(geom_sf(data = analysis_data, aes(color = mass_fdiv_standard, geometry = geometry)), 'Mass FDiv')
std_mass_fdiv_analysis_geo_plot
std_mass_fdiv_analysis_data = model_data(analysis_data[!is.na(analysis_data$mass_fdiv_standard),], 'mass_fdiv_standard')
std_mass_fdiv_analysis_predictors = std_mass_fdiv_analysis_data[,-1]
std_mass_fdiv_analysis_interp = VSURF(x = std_mass_fdiv_analysis_predictors, y = std_mass_fdiv_analysis_data$mass_fdiv_standard)
Thresholding step
Estimated computational time (on one core): 13.2 sec.
|
| | 0%
|
|====== | 5%
|
|============= | 10%
|
|=================== | 15%
|
|========================== | 20%
|
|================================ | 25%
|
|======================================= | 30%
|
|============================================= | 35%
|
|==================================================== | 40%
|
|========================================================== | 45%
|
|================================================================ | 50%
|
|======================================================================= | 55%
|
|============================================================================= | 60%
|
|==================================================================================== | 65%
|
|========================================================================================== | 70%
|
|================================================================================================= | 75%
|
|======================================================================================================= | 80%
|
|============================================================================================================== | 85%
|
|==================================================================================================================== | 90%
|
|=========================================================================================================================== | 95%
|
|=================================================================================================================================| 100%
Interpretation step (on 29 variables)
Estimated computational time (on one core): between 1.2 sec. and 25.8 sec.
|
| | 0%
|
|==== | 3%
|
|========= | 7%
|
|============= | 10%
|
|================== | 14%
|
|====================== | 17%
|
|=========================== | 21%
|
|=============================== | 24%
|
|==================================== | 28%
|
|======================================== | 31%
|
|============================================ | 34%
|
|================================================= | 38%
|
|===================================================== | 41%
|
|========================================================== | 45%
|
|============================================================== | 48%
|
|=================================================================== | 52%
|
|======================================================================= | 55%
|
|============================================================================ | 59%
|
|================================================================================ | 62%
|
|===================================================================================== | 66%
|
|========================================================================================= | 69%
|
|============================================================================================= | 72%
|
|================================================================================================== | 76%
|
|====================================================================================================== | 79%
|
|=========================================================================================================== | 83%
|
|=============================================================================================================== | 86%
|
|==================================================================================================================== | 90%
|
|======================================================================================================================== | 93%
|
|============================================================================================================================= | 97%
|
|=================================================================================================================================| 100%
Prediction step (on 7 variables)
Maximum estimated computational time (on one core): 2.1 sec.
|
| | 0%
|
|================== | 14%
|
|===================================== | 29%
|
|======================================================= | 43%
|
|========================================================================== | 57%
|
|============================================================================================ | 71%
|
|=============================================================================================================== | 86%
|
|=================================================================================================================================| 100%
names(std_mass_fdiv_analysis_predictors[,std_mass_fdiv_analysis_interp$varselect.interp])
[1] "core_realm" "abs_latitude" "city_avg_temp" "city_avg_min_monthly_temp"
[5] "city_avg_max_monthly_temp" "city_avg_ndvi" "city_avg_rainfall"
std_mass_fdiv_analysis_formula = formula_from_vsurp(std_mass_fdiv_analysis_predictors, "mass_fdiv_standard", std_mass_fdiv_analysis_interp)
std_mass_fdiv_analysis_result <- model_average(std_mass_fdiv_analysis_formula, std_mass_fdiv_analysis_data)
Fixed term is "(Intercept)"
std_mass_fdiv_analysis_result_table = model_summary(std_mass_fdiv_analysis_result)
std_mass_fdiv_analysis_result_table
std_mass_fdiv_analysis_pred_plot = plot_vsurp_result(std_mass_fdiv_analysis_result_table)
std_mass_fdiv_analysis_pred_plot
pred_legend <- ggpubr::get_legend(
# create some space to the left of the legend
std_hwi_fdiv_analysis_pred_plot + theme(legend.box.margin = margin(0, 0, 0, 0)) + guides(colour=guide_legend(ncol=2)) + labs(color = "Predictor type")
)
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
geo_legend <- ggpubr::get_legend(
# create some space to the left of the legend
std_mass_fdiv_analysis_geo_plot + theme(legend.box.margin = margin(-80, 0, 0, 12), legend.title.position = "top", legend.key.width = unit(10, 'mm')) + labs(color = "Standardised response")
)
legend = plot_grid(
geo_legend,
pred_legend,
nrow = 1
)
legend
plot_grid(
plot_grid(
std_mntd_analysis_geo_plot + theme(legend.position="none"),
std_mntd_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
nrow = 1
) + draw_label("MNTD", size = 16, angle = 90, x = 0.01, y = 0.5),
plot_grid(
std_gape_fdiv_analysis_geo_plot + theme(legend.position="none"),
std_gape_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
nrow = 1
) + draw_label("Beak Width", size = 16, angle = 90, x = 0.01, y = 0.5),
plot_grid(
std_hwi_fdiv_analysis_geo_plot + theme(legend.position="none"),
std_hwi_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
nrow = 1
) + draw_label("HWI", size = 16, angle = 90, x = 0.01, y = 0.5),
plot_grid(
std_mass_fdiv_analysis_geo_plot + theme(legend.position="none"),
std_mass_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
nrow = 1
) + draw_label("Mass", size = 16, angle = 90, x = 0.01, y = 0.5),
legend,
nrow = 5
)
`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?
ggsave(filename(FIGURES_OUTPUT_DIR, 'process_response.jpg'), width = 3000, height = 3200, units = 'px')
ggplot(analysis_data, aes(x = beak_width_fdiv_standard, y = mntd_standard, colour = core_realm)) +
geom_point() +
ylab("MNTD") +
xlab("Beak Width FDiv") +
theme_bw() + labs(color = "Realm")
ggplot(analysis_data, aes(x = hwi_fdiv_standard, y = mntd_standard, colour = core_realm)) +
geom_point() +
ylab("MNTD") +
xlab("HWI FDiv") +
theme_bw() + labs(color = "Realm")
ggplot(analysis_data, aes(x = hwi_fdiv_standard, y = beak_width_fdiv_standard, colour = core_realm)) +
geom_point() +
ylab("Beak Width FDiv") +
xlab("HWI FDiv") +
theme_bw() + labs(color = "Realm")
mntd_fdiv_analysis = analysis_data %>%
dplyr::select(city_id, mntd_standard, hwi_fdiv_standard, beak_width_fdiv_standard, mass_fdiv_standard) %>%
left_join(community_summary) %>%
mutate(urban_pool_perc = urban_pool_size * 100 / regional_pool_size)
Joining with `by = join_by(city_id)`
mntd_fdiv_analysis
ggpairs(mntd_fdiv_analysis %>% dplyr::select(mntd_standard, hwi_fdiv_standard, beak_width_fdiv_standard, mass_fdiv_standard, regional_pool_size, urban_pool_size, urban_pool_perc), columnLabels = c('MNTD', 'HWI FD', 'Bk FD', 'Mss FD', 'Region Rich.', 'Urban Rich.', '% Urban'))
ggsave(filename(FIGURES_OUTPUT_DIR, 'appendix_standarised_correlation.jpg'))
Saving 7.29 x 4.51 in image